arXiv: 1503.00478vl [cond-mat.stat-mech] 2 Mar 2015 


Memory Effects in Schematic Models of Glasses Subjected to Oscillatory Deformation 


Davide Fiocco, 1 Giuseppe Foffi, 2 ’^] and Srikanth Sastry 3,4, [j] 

1 Institute of Theoretical Physics (ITP), Ecole Polytechnique Federate de Lausanne (EPFL), 1015 Lausanne, Switzerland 

2 Laboratoire de Physique de Solides, UMR 8502, 

Bat. 510, Universite Paris-Sud, 91405 Orsay, France 

3 Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur Campus, Bangalore 560 064, India 

4 TIFR Centre for Interdisciplinary Sciences, 21 Brundavan Colony, Narsingi, 500075 Hyderabad, India 

(Dated: March 3, 2015) 

We consider two schematic models of glasses subjected to oscillatory shear deformation, motivated 
by the observations, in computer simulations of a model glass, of a nonequilibrium transition from 
a localized to a diffusive regime as the shear amplitude is increased, and of persistent memory 
effects in the localized regime. The first of these schematic models is the NK model, a spin model 
with disordered multi-spin interactions previously studied as a model for sheared amorphous solids. 

The second model, a transition matrix model, is an abstract formulation of the manner in which 
occupancy of local energy minima evolves under oscillatory deformation cycles. In both of these 
models, we find a behavior similar to that of an atomic model glass studied earlier. We discuss 
possible further extensions of the approaches outlined. 


I. INTRODUCTION 

Most mechanical systems are constantly subject to de¬ 
formation during their life time. Under certain condi¬ 
tions, such deformations can profoundly alter the micro¬ 
scopic and macroscopic properties such as, in metallurgy, 
the strength (work hardening) or the ductility (strain 
softening). Understanding the effect of such deforma¬ 
tions is clearly of great interest from a practical point 
of view and, at the same time, poses many fundamen¬ 
tal questions [T]. Among different systems, amorphous 
materials are particularly challenging from a conceptual 
point of view, due to the lack of microscopic long-range 
order [2HB|. They are, however, very suitable for many 
practical applications and arise in many context in na¬ 
ture. A prominent example is the case of metallic glasses, 
which are disordered solids that are obtained by fast cool¬ 
ing metallic alloys that have been especially designed to 
avoid crystallization and remain amorphous [7HH]. Other 
widely studied examples, in addition to the most famil¬ 
iar class of molecular and polymeric glass formers, are 
colloidal suspensions maim, foams [12], granular pack¬ 
ings ns ns and biological assemblies such as the cy- 
toskeleton [15 . 

Insight into the behavior of glass formers can be ob¬ 
tained by a computational investigation of their energy 
landscape jTHUTS . The same idea can be applied to me¬ 
chanical deformation of glasses. In this case, the evolu¬ 
tion of the local energy minima, or inherent structures , of 
the model glass, can be followed using a protocol referred 
to as athermal quasi-static (AQS) deformation [20]. It 
has been shown, for example, that systems tend to visit 
energy basins with energies typical of high temperature 
m under shear deformation up to large strains, while 
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under a finite amplitude cycle of back and forth defor¬ 
mation, both rejuvenation (increase in the inherent struc¬ 
ture energy) and over-aging (decrease in the energy) ef¬ 
fects can be observed [221 . 

Recently, it has been shown that very rich and inter¬ 
esting behavior arises when a model glass is subjected 
to repeated cyclic deformations at zero temperature, i.e. 
following the AQS protocol. For a model binary mixture 
glass with particles interacting with the Lennard-Jones 
potential [23] , we have shown that as the amplitude of the 
oscillations increases, the system undergoes a transition 
from a quiescent or localized regime to a diffusive regime. 
In the former, after a short transient, the system remains 
localized in the same energy minimum at the end of each 
cycle, while in the latter state it diffuses in configura¬ 
tion space [21]. This transition from localization to dif¬ 
fusion occurs at a critical amplitude 7 C , as has also been 
reported by other authors on similar systems [25, 26], 
in recent experiments |22 d] as well in other theoret¬ 
ical and computational works HE 133]. As discussed 
in |24| . the transition observed at q c can be identified 
with the yielding transition under steady strain, where 
irreversible behavior sets in. This transition resembles 
strongly the dynamical transition from a reversible to 
an irreversible state that has been found in dilute non- 
Brownian colloidal particle suspensions [31] and granular 
systems H3H3HI33], but the similarities and differences in 
the yielding behavior of these systems merits further in¬ 
vestigation. In that case, however, the reversible states, 
analogous to our localized states, are a consequence of 
the intrinsic reversibility of the low-Reynolds number hy¬ 
drodynamics and they disappear above a certain critical 
threshold when the interactions between the particles sets 
in. 

The existence of reversible or localized states, that re¬ 
main unchanged under the effect of oscillatory deforma¬ 
tion, implies that the system remains indefinitely in a 
given state that corresponds to the amplitude of the os¬ 
cillatory deformation imposed. In principle, this prop- 



2 


erty can be used to encode memory that can be read 
at a later stage. Keim an Nagel have used a compu¬ 
tational model, first introduced in eh. to demonstrate 
that such a memory encoding is possible for a model of 
dilute colloidal particles [32]. In the case of localized 
states, the system retains memory of the value of the 
amplitude that has been used during the training phase 
and this value can be read by analyzing the displacement 
during a deformation cycle, i.e. the reading phase. We 
have shown that surprisingly analogous behavior can be 
observed for a model atomic glass [33], though with im¬ 
portant differences. Glasses possess a complex energy 
landscape and even when the microscopic configuration 
remains invariant at the end of a full cycle of deforma¬ 
tion, the system may move through a complex periodic 
orbit of energy minima during the cycle. In addition, by 
alternating strain cycles of different amplitudes, multiple 
memories of such amplitudes can be encoded. Such a 
possibility exists also in the model colloidal suspension 
but in that case the ability to encode multiple memories 
is transient [34] , and is observed when the system is not 
in the fully trained, stationary, regime, and requires the 
addition of noise to be made persistent for a large number 
of training cycles. Understanding the reasons why these 
(and other) significantly different systems exhibit very 
similar memory effects, and the nature of the universal¬ 
ities and differences is an interesting subject for further 
investigations. In this paper, we describe results concern¬ 
ing two simple, schematic models, which have been inves¬ 
tigated to elucidate the characteristics and mechanisms 
of the non-equilibrium transition and memory effects de¬ 
scribed above. These simplified models are investigated 
in order to understand what are the fundamental mech¬ 
anisms that can originate the phenomenology that has 
been observed for the LJ system. The first is a lattice 
model, the NK model, that we have already briefly dis¬ 
cussed earlier in the context of memory effects [24] • The 
NK model is a spin model with multi-spin interactions 
that is designed to generate energy landscapes with tun¬ 
able roughness, and also to mimic the effects of impos¬ 
ing shear deformation. The second model, the transition 
matrix (TM) model, aims to capture key aspects of the 
changes induced by cyclic shear deformation by specify¬ 
ing rules to construct a transition matrix that maps the 
set of local energy minima onto itself as a result of a sin¬ 
gle cycle. Studying the evolution of the set of occupied 
local energy minima as a function of the number of cy¬ 
cles allows us, in principle, to understand what are the 
key features of the transition matrix that determines the 
observed behavior. 


The paper is organized as follows. In fection II| and [sen] 
Ition IIll we introduce the two models that are the main 
subject of this paper. In section IV we discuss the re¬ 
sults that concern the existence of a dynamical transition 
while in |section V| we discuss memory effects. Finally, we 
summarize our results and conclusions in lsection VII 


II. THE NK MODEL 

The NK model is defined on a set of spins interacting 
according to a Hamiltonian that incorporates a parame¬ 
ter 7 that mimics the effect of shear strain. It was studied 
by Isner and Lacks [36] as an analog system to the model 
atomic glass former in which overaging and rejuvenation 
were observed | 22 j under a single deformation semicycle. 

We consider N (even) lattice sites occupied by spins 
to,; that can take either the values 0 or 1 . 

{toi,to 2 , ... ,to, ; , ... ,m N } £ { 0 , 1 }^ ( 1 ) 

Furthermore, in order to prevent the system getting 
trapped around low energy configurations with aligned 
spins, we limit the space of allowed configurations to 
those that satisfy the constraint 

= y ( 2 ) 

i 

(this is equivalent to taking the set of states of constant 
magnetization 0 in the Ising model). There are ( N / 2 ) 
such configurations, and we define two NK configurations 
as adjacent if one is turned into the other by swapping the 
values at two sites i and j such that to, ^ mj. Since one 
may choose one of N/2 spins with to = 1 for swapping 
with any one of N/2 spins with m = 0, each configuration 
has N 2 / 4 configurations that are adjacent to it. 

We introduce: 

• An ordered list of K “neighbors” for each *-th spin, 
specified by the map J: 

mi -4 } (3) 

The choice of the list of neighbors for a given spin 
is random. 

• Two functions a and b connecting the set {0, 1} K+1 
(the set of all the (A'+ 1 (-tuples formed by ones and 
zeros) to the intervals [— 1 , 1 ] and [ 0 , 1 ] respectively 

{0, 1} K+1 A [—1,1] {0, 1} K+1 \ [0,1] (4) 

The correspondence between a given (K + l)-tuple 
and the numerical value is chosen randomly with a 
uniform probability within the respective intervals. 

The energy of the system depends on the values of the 
spins, E = A(toi, ..., uin) according to 

1 N 

E = --^2 [i + sin(27r(a, + 7 &i))] (5) 

z *=l 

where 7 (which we will name, with an abuse of language, 
“shear strain”) is a parameter that can be varied contin¬ 
uously. 

The NK model is known to possess a discrete energy 
landscape whose roughness (a measure of which is given 
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by the number of local minima in the landscape) is tuned 
by the value of the parameter K. The landscape is 
smooth for K = 0 and the roughness is increased as K is 
increased. To see this consider the case 7 = 0 for simplic¬ 
ity. If K = 0, when performing the sum in |Equation 5| one 
sums contributions that are simply 1 + sin (27m,), where 
di can assume only two values in the interval [— 1 , 1 ], de¬ 
pending solely on whether rrii = 1 or 0 , as no spin has 
neighbors. It is clear that the energy increases or de¬ 
creases monotonically with the “magnetization”, [U ,, nii. 
As in our case such quantity is fixed, all the allowed con¬ 
figurations have the same E. This implies that by moving 
from a given configuration to any adjacent one the en¬ 
ergy can’t change and the landscape is thus flat. In the 
case K = N — 1 each spin is a neighbor of any other. 
To estimate the roughness of the landscape, we compare 
the energies of two adjacent structures, that differ by the 
swapping of two spins. Each term in |Equation 5| is af¬ 
fected by such a swap, and so all the ad s that contribute 
to the energy of the two configurations differ. Since the 
ad s are random and uncorrelated, there won’t be correla¬ 
tion between the energies of two adjacent configurations. 
Consequently, the structure of the overall landscape will 
be rough. 

While I\ is able to tune the roughness of the landscape, 
the parameter 7 is able to change continuously the values 
of the energy of the configurations (and thus the overall 
landscape). This is similar to what happens in a glass, 
in which a macroscopic strain is externally imposed. 

The broad features of the NK model are thus similar 
to the atomic glass formers mentioned earlier, with the 
significant difference that the configuration space of the 
NK model is discrete. We next describe how procedures 
such as energy minimization and AQS are performed in 
the NK model. Starting from an initial configuration, 
energy can be minimized by a steepest descent (SD) pro¬ 
cedure. A SD in the NK energy landscape consists of 
moving from a configuration to the adjacent one with 
the lowest E, and iterating this procedure until when no 
move to an adjacent configuration results in a decrease 
in E. Using such a protocol, any configuration can be 
mapped onto a local minimum, i.e. an inherent struc¬ 
ture of the NK landscape. This fact, together with the 
dependence on 7 of the landscape, allows us to define an 
athermal quasi static “deformation” procedure on NK 
systems too: 

1. Take an inherent structure of the NK landscape. 
This can be obtained by any configuration by 
means of the SD algorithm. 

2. Increment the value of 7 in |Equation 5| by a small 
amount dj. The value of d'y should be small enough 
so that the modification of the landscape is slow 
enough and no displacements to adjacent inherent 
structures are “missed” by the AQS dynamics. 

3. Apply the SD procedure to the configuration. 


As AQS can be applied to the NK model, it is inter¬ 
esting to check whether the same phenomenology seen in 
AQS deformation of model atomic glass formers can be 
observed in it. However, one should bear in mind the 
following differences between the NK model and atomic 
glass formers: 

• Since the NK model has a discrete configuration 
space, the AQS dynamics of NK systems is some¬ 
what different to that of atomic systems under the 
AQS protocol. In the NK case, the system occupies 
one given point of the available configuration space 
as 7 is changed, and stays there until it “jumps” to 
another inherent stucture configuration as soon as 
the initial configuration is not an inherent structure 
anymore. In the atomic systems, instead, the con¬ 
figuration of the system continuously changes as 7 
is varied. 

• Due to the discrete nature of the landscape, min¬ 
imization is trickier in the NK case. While for 
atomic glass formers local quantities (e.g. the cal¬ 
culation of a potential energy gradient) allow the 
determination of directions to follow to reach an 
inherent structure, in the NK model all the ener¬ 
gies of adjacent configurations need to be calculated 
in order to choose the adjacent configuration with 
the lowest energy (if there is one). This opera¬ 
tion requires 0(N 2 ) energy calculations to be per¬ 
formed (as the number of adjacent configurations 
scales with N 2 , see above) and it is thus computa¬ 
tionally infeasible for large values of N. 

• While for atomic glass formers, the energy land¬ 
scapes depend on the few parameters defining the 
interaction potential (e. g. the values of the e and 
a parameters for the Lennard-Jones potential) and 
the boundary conditions of the simulation volume, 
the definition of the landscape in the NK case re¬ 
quires the introduction of a much larger number of 
parameters. This is because the lists of neighbors 
specified by J can be realized in many different 
ways and a and b in |Equation 3| and [Equation 4| 
require 2 K values each to be defined. 

Having noted this, we can study the NK model under 
athermal quasi-static deformation, and compare the re¬ 
sults with what is obtained with the model atomic glass 
former studied in |24j (the Kob-Andersen binary mix¬ 
ture with Lennard-Jones interactions [23] , or LJ). Such 
comparison is meaningful, because the two systems share 
important features. At the same time it is not trivial, 
becuase, as explained above, the two models are “suffi¬ 
ciently different” that common qualitative behavior can 
not be taken for granted a priori. 

Before performing simulations of “deformation” on NK 
samples one needs initial configurations. As in the case of 
the LJ systems discussed in [21], one would like to start 
from sets of configurations that differ for some feature 
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FIG. 1. Average energy per site of the inherent structures 
obtained by quenching NK samples (N = 20, K = 10) equili¬ 
brated at different T. 


(like the average potential energy). In this way one can 
check how deformation affects the samples, and how it is 
capable of making them evolve in such a way to forget 
their initial state. The obvious choice is to choose inher¬ 
ent structures corresponding to different temperatures T, 
as we did in our previous study of the LJ case. We thus 
obtain several equilibrated NK configurations at different 
T. We do this for several realizations of the couplings a, 
b and J by means of Monte Carlo sampling and find the 
inherent structures at each temperature T by means of 
steepest descent (SD) minimization. The average energy 
of such inherent structures is plotted in |Figure 1[ showing 
a behavior similar to that observed in j 2 S] for a model 
atomic glass former. We choose temperatures (measured 
in units of ks = 1) T = 0.6 and 1.0 for further anal¬ 
ysis, as they are different enough to be distinguished. 
These temperatures are above the glass temperature of 
the models which is estimated to be around 0.45 {SHI ■ 


III. TRANSITION MATRIX MODEL 

Both the LJ and NK models possess a rugged land¬ 
scape that is modified during the deformation. The be¬ 
havior under oscillatory deformation, in particular, thus 
depends on the detailed features of the landscape. Would 
it be possible to predict the behavior of such models 
qualitatively, without encoding in detail the features of 
the energy landscape, but using a sort of “high-level”, 
abstract description of its evolution? This is what the 
“transition matrix” method (TM) aims to do. 

The starting point of the TM approach is to look at a 
cycle of AQS deformation from a mere “mathematical” 
point of view. In that perspective, an AQS cycle is a cor¬ 
respondence between the set of M inherent structures of 
the energy landscape at 7 = 0 into itself. This is because 


a valid starting configuration is an inherent structure of 
the 7 = 0 landscape, and it is transformed into another 
inherent structure of the same landscape at the end of 
the deformation cycle. Each of these inherent structures 
can be identified by an index i. and associated to a M- 
dimensional vector Rj whose components are all zero but 
for the *-th one, which is set equal to 1. Then, they can 
be taken as starting points of a deformation experiment 
where a single AQS cycle is performed and the inherent 
structure reached at the end is recorded. 

This allows to define a transition matrix P, such that 

PR; = R f ( 6 ) 


where Ri and Rf are the vectors associated to the ini¬ 
tial and final inherent states. Here we list some of the 
properties of the P matrix, which derive directly from 
the features of AQS dynamics: 

• P encodes the entire information about the evolu¬ 
tion of any inherent structure under cyclic defor¬ 
mation, as by using [Equation 6 ] one can determine 
the final inherent structure Rf given any intial one 

R ; . 


• P is a sparse M x M matrix, and P l; j = 1 if and only 
if the state associated to Rj is mapped onto the 
inherent structure Ri in the AQS cycle. The con¬ 
sequence is that all the columns have exactly one 
non-zero entry (equal to one) because each Rj con¬ 
figuration is sent to some other Ri inherent struc¬ 
ture by the AQS cycle. 

• P depends on the value of 7 max, be. on the am¬ 
plitude of the deformation. For very small ampli¬ 
tudes, a sizable fraction of the inherent structures 
will be unchanged under the deformation, because 
an AQS cycle will not be effective at destabilizing 
the starting inherent structures. In this case a given 
structure will often map onto itself through P and 
thus P will be very close to the diagonal unit ma¬ 
trix. In general this won’t be any longer true for 
higher values of 7 max . 


• The determinant of P is, in general, zero. In gen¬ 
eral, more than one structure will map onto the 
same final state Rf, so that in that case P will de¬ 
fine a non-injective function. As the domain of P is 
also its codomain, P is not surjective, so that there 
will be structures that are not arrival configurations 


for any inherent structure (this is illustrated in Fig¬ 


ure 2). This has the consequence that some rows of 


the P matrix, in general, are identically zero, and 
so is the value of det P. 


• M, in general, is a large number. The size of P is 
equal to the number of inherent structures of the 
landscape, and this number, for typical LJ and NK 
landscapes, is large (exponential in N in the first 
case, and equal to (^ 2 ) the secon d)- 
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FIG. 2. Graph representation of an AQS cycle, that has the 
effect of mapping the set of inherent structures of the 7 = 0 
landscape onto itself. 


• The result of a deformation experiment where L 
cycles (rather than just one) are applied to a given 
starting configuration Ri is obtained by applying 
[Equation 6 | repeatedly: 

P l R ; = R f where P L = P • . . . • P (7) 

L times 

A. Classification of states by their transformation 
properties 

A given configuration Rj can transform under the ef¬ 
fect of P in different ways: 

1. PRi = Ri: in this case P has no effect on Ri, 
so that Ri is an eigenvector of P relative to the 
eigenvalue 1 , and Pa = 1 . We name such a Ri an 
absorbing state. 

2. P L Ri = Ri for some L > 1: in this case the oscil¬ 
latory deformation starting from Ri makes it cycle 
through a sequence of states, and after L cycles Ri 
is reached again. Ri an eigenvector of P L relative 
to the eigenvalue 1. We name such a Ri a recurring 
state. 

3. P J Ri = R a bs for some J > 1 , where R a bs is an 
absorbing state. We name such a Ri as mapping to 
absorbing state. 

4. P J R ; = R roc for some J > 1, where R rC c is a 
recurring state. We name such a Ri as mapping to 
recurring state. 

Note that every configuration falls in one of the cate¬ 
gories enumerated above. This can be demostrated as fol¬ 
lows: suppose that a configuration X exists that doesn’t 
fall in any of the categories listed above. If P is applied 


to it, then some other configuration Yi 7 ^ X is obtained, 
else X would belong to the category in |item 1 Repeated 
application of P yields always new states Y 2 , Y 3 ,..., so 
that no absorbing nor recurring states are encountered, 


else X would belong to one of the categories in item 2 


|item 3| or |item 4| After M applications of P, a sequence 
of M + 1 distinct inherent configurations has been gen¬ 
erated, but only M distinct inherent structures exist! So 
the initial hypothesis (the very existence of X) can’t be 
true. We note that in principle, recurring states with 
periods that are of 0 ( 1 ) are qualitatively different from 
and should be distinguished from those with are of 0(M), 
although we do not attempt such analysis here. 


What is the correspondence between this classification 
of states and the absorbing and the diffusive states en¬ 
countered in [M|? Absorbing states of the LJ model 
clearly correspond to the absorbing states of the TM 
model. Diffusing states in the LJ model correspond to 
a subset of the recurring states of the TM picture with 
periods which are large, or more precisely, exponential 
in the number of particles in the system. In fact, even 
though a diffusing LJ system does not seem to revisit 
the same inherent structure as it travels in configura¬ 
tion space, after a large number of oscillation cycles it 
has to, as the number of possible inherent structures is 
finite. Thus the diffusing states of can be viewed as re¬ 
curring states, which just take a very large number of 
cycles to come back to their starting state. Neverthe¬ 
less, a more sophisticated analysis than we attempt here 
should distinguish between diffusing states and recurring 
states with 0 ( 1 ) period cycles. 


B. Construction of the P matrix 


[Equation 6 | tells that the P matrix contains the en¬ 
tire information about the outcome of oscillatory AQS 
deformation, but how can one construct it? Comput¬ 
ing it for LJ or NK systems can in principle be done 
by brute force: one needs to have a list of the inherent 
structures, use each of them as a starting configuration 
for a shear deformation cycle and see which structures 
they eventually reach at the end of the cycle. This idea 
is easier to apply for the NK model than in the LJ case: 
in the former one can (at least in principle) enumerate 
all the ( N / 2 ) allowed configurations, minimize each and 
every one of them so to get all the inherent structures 
of the landscape; in the latter, which has a continuous 
energy landscape and a continuous set of configurations, 
the determination of the local minima of the landscape 
is a less trivial task[3S|. The infeasibility of a brute force 
approach is the reason why one would like to construct 
P by less expensive means, albeit in an approximate or 
schematic way. To do so, we make a series of observa¬ 
tions and assumptions about the evolution of the energy 
landscape, that allow to construct P. 
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1. Assumptions for constructing transition matrices 

The transition matrices that we construct are based on 
some assumptions on the dynamics of inherent structures 
in the course of a deformation cycle: 

1. During the deformation (at non-zero values of 7 ) 
the number of energy minima is assumed to be al¬ 
ways equal to M, no matter the value of 7 . The 
number of local minima present in the landscapes 
of LJ and NK models will, in general, weakly[40] 
depend on the value of 7 . 


0.750 

-0.755 


• • T = 0.6, - 0.3 

T = 0.6, 7™* = 0.4 

• • T = 0.6, 7 max - 0.5 

-0.760 


. T = 0.6, 7 max - 0.6 

• • T = 0.6, 7 max = 0.8 

• • T = 0.6, 7 max - 1.0 
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T = 1.0, Jmax - 0.4 
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a a T — 1.0, 'Ymax - 0.8 

-0.775 


* T = 1.0, Jmax = 1.0 

* * A V 

• ..... 


O .785 

y 

0 — 
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7< 


2. Minima are destabilized by changing 7 , i.e. that 
some of them are destroyed by the deformation. 
As the number of structures is assumed to be con¬ 
served (see above) to each inherent structure de¬ 
struction corresponds the creation of a new one. 

3. The probability per unit strain of an inherent struc¬ 
ture to be destabilized is assumed to be the same 
for all the M structures and equal to a value r in¬ 
dependent of 7 . 

4. A system resides in a given inherent structure un¬ 
til such inherent structure is destabilized. When 
this happens, the system jumps to another inher¬ 
ent structure of the deformed landscape. For sim¬ 
plicity such a structure is assumed to be picked at 
random in the landscape (whereas, in a realization 
of the NK or LJ models, a system will land on a 
structure which is not far from the starting one in 
the space of configurations). 

In addition, the model relies on two facts that are true 
in general: 

5. As a deformation semicycle brings the system from 
0 up to 7 ma x and then back to the undeformed land¬ 
scape, there is a symmetry in the structures that 
are created and destroyed as 7 is incremented from 
0 to 7 max and those that are created/destroyed as 
7 is reduced back to 0 in the second part of the 
semicycle. In fact, if a structure R is destroyed 
when incrementing 7 above some value 7 *, the same 
structure R will be created at 7 * as the deforma¬ 
tion is reversed. The converse is true for a structure 
S that is created in the first half of the semicycle. 

6 . The matrix P can be viewed as the product of P + 
and P -, the matrices that describe the two semicy¬ 
cles (one denoted by positive, the other by negative 
strain 7 ) that form a full oscillation cycle. 

A description of how these assumptions and observa¬ 
tions are combined to calculate an approximation of the 
P matrix is presented in Appendix |A| 


FIG. 3. Potential energy per site as a function of 7 acc , for 
different initial effective temperatures and different deforma¬ 
tion amplitudes 7 maa ,, setting N = 40, K = 10. Data refer 
to configurations with 7 = 0 . For large values of 'y ma x the 
energy fluctuates around some value which depends on 7 max 
only. At small 7 mox , instead, the plateau value of the energy 
depends on the effective T of the initial configuration. In this 
respect what is observed here qualitatively resembles what 
has been found in the atomistic glass model in 124) . 


IV. RESULTS: DYNAMICAL TRANSITION 

A. Dynamical transition under oscillatory 
deformation: the NK model 

We consider N = 20,40, 80 and K = 10, and ss 200 
instances of the couplings and obtain 3-4 equilibrated 
configurations at T = 0.6 and 1.0 for each of such in¬ 
stances. The corresponding inherent structures obtained 
by SD are deformed by increasing 7 in |Equation 5| in 
steps of (I 7 = 0.005 (for all TV) and performing SD at 
each step. The parameter 7 is varied in the interval 
hlmoii Ttnoi] in a triangle wave fashion, exactly as 
in the AQS simulations seen in Ref. [Ml and (33 The 
values of E and the configuration are recorded whenever 
7 = 0, i.e. at intervals of 2"/ max . Plots of £ as a 
function of the accumulated strain 7 acc are shown in 
|Figure 3| for the TV = 40 case. Similarly to what is 
observed in the LJ case, for small values of 7 max the 
energy reaches a plateau which depends on both 7 max 
and the initial effective T. For higher values of 7 max , all 
samples, regardless their initial effective T, reach plateau 
depending on 7 max only. In this respect, the NK model 
is able to reproduce qualitatively the same behavior 
found in LJ systems [Hj, where samples forget about 
their initial preparation if the oscillation amplitude 
exceeds some value j c . 

Changes in configurations under deformation can 
be studied for NK configuration as well, by looking at 
the distance between configurations before and after the 
application of deformation cycles. Consider two configu¬ 
rations, R( 7 a CC ) and R( 7 acc )> obtained for values of the 
accumulated strain equal to 7 acc and 7 acc respectively, 
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with "f acc < "face- Their distance can be expressed using 
the Hamming definition [4T] d: 


diriacc - lacc) = 


Coi + Cio 

N : 


( 8 ) 


where coi (cio) is the number of occurrences such that 
the i-th component of R( 7 QC c) and the *-th component of 
R(7 ace) are respectively equal to 0 and 1 (to 1 and 0). We 
thus pick a configuration r 1; choosing a large enough 7 acc 
so that the corresponding E in |Figure~3] has relaxed to 
a steady state. We then compute the Hamming distance 
from it for configurations reached for increasing values of 
7 occ , and plot it in |Figure~4| 

The average Hamming distance measured starting 
from a reference configuration in the steady state quickly 
reaches a constant value for increasing 7 acc , rather than 
show a linear increase with time as the mean squared 
displacement does for the atomic glass former [21] . How¬ 
ever, the Hamming distance, like the overlap [12] or dis¬ 
similarity [43j used in the context of atomic liquids, is 
a bounded quantity and should be treated as a correla¬ 
tion function rather than a measure of mobility. In the 
present case, given the constraint on the total magneti¬ 
zation, the Hamming distance has an average value of 
0.5 for uncorrelated configurations. Thus, the value of 
the Hamming distance itself serves as a metric of local¬ 
ization of the configurations under shear. 

By plotting the average Hamming distance as a func¬ 
tion of 7 max for different values of N (see Figure 4), 


one observes that the capability of the system to move 
away from the reference configuration increases sharply 
at some value 7 C , which is roughly N independent. The 
sharpness of the transition, moreover, is seen to increase 
with N. These data thus seem to confirm that a tran¬ 
sition at some oscillation amplitude y c from a localized 
regime to a diffusive one exists in the NK model, simi¬ 
larly to what is observed in our LJ Systems in Ref. |2D 
and models of particle suspensions studied in [54] . 




FIG. 4. (top) Hamming distance as a function of the accu¬ 
mulated strain measured from reference configurations with 
7 a cc > 7 acc, where 7 acc marks the reaching of the plateau of 
the energy in |Figure 3| IV = 40 and K = 10. The behavior is 
not diffusive, but the higher the 7 max, the further the systems 
are able to move away from the reference configuration. The 
Hamming distance can be modeled with a constant function 
of 'Yacc — "face- (bottom) Average value of the Hamming dis¬ 
tance as determined by averaging data like that in |Figurc 4| 
for different N and setting K = 10. The ability of the sys¬ 
tem to diffuse away from a reference configuration increases 
strongly with 7 max at 7 max « 0.5. The transition becomes 
sharper as N is increased. 


B. Dynamical transition under oscillatory 
deformation: the TM model 


Since the transition matrix model does not contain 
any real or configuration space distance information, a 
suitable measure of localization has to chosen for this 
model based on information about the fate of individual 
inherent structures under repeated operation of the 
transition matrix. The information contained in P can 
be used to distinguish states that are absorbing or map¬ 
ping to absorbing states from those that are recurring 
or mapping to recurring ones (see the definitions given 
in section III A). This information, in turn, can be used 


to gather information about the dynamics under AQS 
deformation: if absorbing states dominate, systems are 
likely to be trapped into them, similarly to what is 
observed in the LJ and NK models below 7 C ; if recurring 
states dominate, systems have the capability of exploring 


the configuration space before returning to the same 
state point, analogously to what happens in the other 
models above y c . 

The naive way to determine how inherent struc¬ 
tures behave under oscillatory deformation is to apply P 
repeatedly to each of them. After some applications of 
P, the structures will transform into states of the kind 
R a bs such that PR a b s = R a bs (absorbing), or states 
of the kind R rec such that P L R rec = Rrec (recurring). 
This procedure is however too expensive, because one 
needs to calculate the trajectory of each and every of 
the M structures by means of matrix multiplication. 
A way to overcome this is treating P as an adjacency 
matrix, and constr ucting the directed graph[44j G 
associated to it (see Figure 5). G will be a directed 
graph whose outdegree is 1 , as each structure maps onto 
one and only one configuration through P. In general, 































G will possess several connected components^;. Each 
of these either contains a self-loop or not. Connected 
components containing a self-loop are those that contain 
an absorbing state, which is a node that is connected to 
itself via the self-loop. All the other nodes are connected 
to it, and thus represent states mapping to the absorbing 
state. Connected components not containing self-loops 
must contain a loop, and their nodes thus represent 
recurring states or states mapping to recurring states. 
By examining the graphs one can count the number R 
of recurring (and mapping to recurring) states simply 
by counting the number of nodes of the connected 
components of G which do not possess self-loops. The 
number of absorbing (and mapping to absorbing) states 
will be given by M — R. 


We obtain P with the procedure described in section A 


using Python and the support for sparse matrices within 
the library SciPy 03 ■ We then extract G and its con¬ 
nected components using NetworkX [35] . The connected 
components associated to recurring states can be easily 
filtered because they don’t contain self-loops. We do so 
for matrices P with M = 10 4 , 10 5 and 10 6 , setting (ten¬ 
tatively) the probability for an inherent structure to be 
destabilized per unit strain to r = 0.04 and plot the aver¬ 
age fraction of recurring (or mapping to recurring) states 
as a function of the 7 max averaging on ss 800, 200, 50 


matrices respectively. The result is shown in Figure 6 


The curves in |Figure 6 | can be modeled by the fitting 
function: 



Data in |Figure~ 6 | and the form of |Equation 9| show that 
the TM model shows a sharp increase in the number of 
states mapping to non-absorbing states as the “oscilla¬ 
tion amplitude” is increased beyond some value 7 c , sim¬ 
ilarly to what has been observed in the case of LJ mix¬ 
tures and of the NK model above. For this reason, one 
can reasonably believe that while crude, the TM model 
capture enough details to describe qualitatively the tran¬ 
sition from a “localized” regime (where absorbing states 
prevail) to a “diffusive” one (where recurring states dom¬ 
inate) observed in particle models. Moreover, the transi¬ 
tion is observed to be sharper for higher values of M, with 
the parameter a in |Equation 9| increasing with increasing 
M. Opposite to what is observed in LJ systems, how¬ 
ever, the value of 7 c is seen to increase with the system 
“size” M. It would be useful to address the question of 
whether a sharp transition arises in the thermodynamic 
limit rigorously [4SJ- 


V. RESULTS: MEMORY EFFECTS 

In earlier work [35] the memory effects in a model 
atomic glass former (LJ) and the NK model were studied, 
and shown to be similar, and we showed that both sys- 




FIG. 5. (Top) The output of the TM model is a map (which 
depends on 7 max ) of the set of inherent structures onto itself. 
The same map can be interpreted as an adjacency matrix for 
a directed graph. (Bottom) The resulting graph is a collection 
of 1-trees [46]. Each of these 1-trees can contain a self-loop or 
not. If it does, then the vertex with the self-loop represents 
an absorbing state, and all the vertices in the 1-tree to which 
it belongs represent states mapping to that absorbing state 
via AQS dynamics. If the 1-tree does not contain self-loops, 
then its vertices are associated either to recurring states or 
states that map onto recurring states in the AQS dynamics. 
By counting the number of vertices in the two kinds of 1-trees 
(those with self-loops and those without), one can thus deter¬ 
mine the fraction of states that are absorbing (and mapping 
to absorbing) or recurring (and mapping to recurring). 


terns were capable of encoding persistent multiple mem¬ 
ories. We thus focus here on the memory effects found 
in the TM model, but compare them with the behav¬ 
ior for the cases studied earlier. The procedure to probe 
memory effects in the TM model is fairly different with 
respect to that followed in the LJ and NK cases. First of 
all, we generate P matrices for different values of 7 max , 
following the procedure described in |section A] Once a 
71 is chosen, we consider the configurations trained by 
N cyc oscillations (equivalent to 7 acc = 47 !lV cyc ). These 
are configurations for which the corresponding rows of 
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FIG. 6 . Overall fraction of recurring states and mapping to 
recurring states as a function of 7 max for different values of 
M, obtained by analysis of the graphs associated to the tran¬ 
sition matrices P generated by the TM model. The number 
of recurring states increases strongly beyond some value 7 C 
which increases as M increases. The data can be described 
fairly well with the model in Equation 9 (black curves), with 


a sharpness (dictated by the parameter a) which increases for 
larger M. 


the matrix P^[ avc have at least one non-zero entry, where 
P 7l is the matrix associated to the deformation up to 71 . 
The configurations that, after the N cyc of training, cor¬ 
respond to rows with all zero entries in the deformation 
matrix P^[ ovc cannot be reached by any further transfor¬ 
mation and cannot be considered. The remaining ones 
are those that can be still trasformed under the effect 
of the deformation. We call to the set of such states 
An c . To probe the behavior of such states under a sin¬ 
gle reading cycle of amplitude j r , we check whether they 
are absorbing cycles for 7 r , i.e. we verify if the condition 
P 7i R = R is satisfied, for each R state in An c c . The 
validity of this condition is easy to check, because it holds 
true if and only if the matrix element on the diagonal of 


P 7r . corresponding to such states is equal to one. In Fig¬ 


ure 7 we plot the fraction of non-absorbing states (i.e. R 


states in Am that don’t meet the condition P. R = R) 
for different values of 7 r , obtained by studying the states 
belonging to a space with M = 10000 structures with 
probability r = 0.04 to be destabilized and “trained” by 
a different number of cycles of amplitude 71 =60. 

The results shown in |Figurc 7| indicate that the TM 
model displays single memory effects very similar to the 
LJ and NK models. Strictly speaking, the plot in |Fig-| 
ure 7[bottom) is not equivalent to those in Figure 7~| (top) 
and Figure 7| center): in the LJ and NK models a no¬ 
tion of distance exists between configurations, so that 
one is able to quantify the displacement experienced by 
the samples during the reading phase (using the MSD 
and the Hamming distance respectively); such informa¬ 
tion is not available for the TM model, which however 
offers very similar information: we can answer whether 


states are left unchanged or are modified by a reading 
cycle of amplitude 7 ,. [50]. 


A. Multiple memories 

The training can be modified so that it consists of al¬ 
ternated repetition of cycles of amplitude 71 and 72 , like 
in Figure 8 | The rationale behind this protocol is to 


encode multiple memories in our samples, and be able 
to read the values of the different training amplitudes 
in the reading phase. The cycles in 7 have the form 
0 —» 71 —t 0 —> —71 —l 0 — l 72 —1 0 — > —72 —► 0. We 

repeat such cycle N cyc times, so that after the training 
the sample as been subjected to an accumulated strain 
lacc = 4(71 + "f 2 )N C y C - This can be straightforwardly 
done in the LJ and NK cases, whereas (as described above 
in the case of single memory) a different scheme must be 
adopted with the TM model. We compare results for the 
LJ, NK and TM models. 

For the LJ model, we choose 71 = 0.06 and 72 = 0.04 
and train samples of the same size and initial effective 
temperature as those trained with a single amplitude by 
performing N cyc on them. We then take copies of the 
trained samples and subject them to a reading cycle of 
amplitude 7 r . As above we measure the MSD of the con¬ 
figurations as a function of j r . As it can be seen from 
Figure 9 [ top) , the MSD has two kinks in correspondence 
of 71 and 72, which are both visible for sufficiently high 
N cyc . In addition, for a high number of N cyc , the MSD 
curve converges to a curve showing clearly the trace of 
the two training amplitudes. By looking at the data, it 
is reasonable to assume that this will be true for an ar¬ 
bitrarily large number of N cyc , so that the two memories 
will be persistent for 7 acc —>• 00. 

For the NK model, we choose 71 = 0.06 and 72 = 
0.04 and train the samples for different N cyc with the 
same N, A', initial effective temperature and values of 
the couplings of those trained with a single amplitude. 
Again, we use the Hamming distance d as a measure of 
the change in configurations in the reading phase and 
plot d for different values of j r in Figure 9 ] ( center). As 
in the LJ case, the training amplitudes can be read by 
looking at kinks (discontinuities in the first derivative) 
in the plot, and again they appear persistent in the limit 
7acc ^ OO- 

In the case of the TM model the “trained” configu¬ 
rations are those with the same index of the non-empty 
rows of the matrix (P 72 P 7l ) ,v ' :,Jc , where P 7l and P 72 are 
the matrices associated to the deformation up to 71 and 
72. We refer to the set of such states as -Bjv c „ c . Ex¬ 
actly as in the reading of single memories, we compute 
the fraction of non-absorbing states (i.e. states in B^ c c 
that don’t meet the condition P lr R = R) for different 
values of ”f r ) as a function of 'fmax- From the analysis 
of Figure 9 (bottom), it is clear that a double memory 
can be encoded in an ensemble of structures in the TM 
model. 
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VI. CONCLUSIONS 

In this paper we have investigated the role of oscilla¬ 
tory deformation in two toy models that have been intro¬ 
duced to mimic the behavior of a model of glass former 
that we studied before {24[ [35]. In particular we were 
interested in the presence of a dynamical transition from 
a localized to a diffusive state upon increasing the ampli¬ 
tude of the oscillations and in the encoding of memory in 
the system when a proper training protocol is followed. 
The first model that we discussed is the NK model, a 
disordered spin model that has been introduced to inves¬ 
tigate the deformation in glasses [55J ■ The main property 
of the NK model is an energy landscape with tunable 
roughness. We have shown in our previous work that 
in this model we could induce a memory effect by oscil¬ 
latory training cycles in a similar way to what can be 
achieved in LJ glass formers [35] and colloids [34.;. Here 
we have discussed in more detail the nature of this en¬ 
ergy landscape and we have shown that under oscillatory 
deformations, the NK model presents a dynamical tran¬ 
sition at a critical amplitude of strain that is very similar 
to the one that is found in glass formers [23]. 

We also introduced a second model, the TM model. This 
model is based on a matrix approach that represents a 
further abstraction with respect to the NK model. The 
main idea is that, after a deformation cycle, when the 
sample returns to a state of zero deformation, the po¬ 
tential energy landscape is fixed. Consequently, one can 
map the act of deforming into the action of a transition 
matrix that changes the occupied minima of the land¬ 
scape that are represented in a vector space. We have 
shown that, by imposing a few hypotheses, a reasonable 
transition matrix can be built. Despite its simplicity, we 
have shown that this model is capable of reproducing 
both the dynamical transition and the memory effects 
that we have found in the LJ and NK models [24] [35] . 

In conclusion, the NK and TM models are interesting 
for two main reasons. Firstly, they clarify what are the 
essential model features behind the dynamical transition 
observed in the LJ model. The NK model, for example, 
proves that a continuous configuration space and energy 
landscape is not necessarily required. The TM model 
goes even further and shows that the observed phenom¬ 
ena occur when the transition matrix is constructed fol¬ 
lowing a simple set of rules. Having found that such a 
set of basic ingredients exist, one can imagine that anal¬ 
ogous dynamic and memory effects could be observed in 
a wide class of systems, such as spin models or soft mat¬ 
ter materials. Secondly, these simplified models could be 
the base of more fundamental approaches to treat the 
athermal deformation of physical system, in which the 
dynamics is deterministic. The TM model, for example, 
could lay down the basis for a new framework to model 
oscillatory deformation in materials. 

The two models that we have described here make some 
fundamental simplifications on the properties of the PEL. 
The TM model assumes that the states visited during 


deformation are completely uncorrelated (as there is no 
notion of space or distance in the TM model), while the 
NK introduces a PEL of controlled roughness but with 
no structural information, as it is generated from a se¬ 
ries of random couplings. This is fundamentally differ¬ 
ent from the full KA model that presents a PEL that is 
the result of a two-body interaction potential between its 
constituents. A present limitation of the TM model, how¬ 
ever, is that it doesn’t incorporate any organization be¬ 
tween the different inherent states. This is clearly not the 
case for the LJ systems, whose PEL has a complex struc¬ 
ture and spatially correlated particle motion between ad¬ 
jacent inherent structures is expected under deformation. 
In the future, it would be interesting to see if some degree 
of correlation can be included in this model. Another in¬ 
teresting question for future work is the identification of 
precise conditions for persistent memory as seen in the 
NK and TM models, as opposed to transient memory 
seen in other model systems [33] . 
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Appendix A: Details of the construction of P can see that 


P = P_P + 


(Al) 


Here we describe how the assumptions listed in |sec- 


tion IIIB 1 can be used to construct P. First of all, one 


so that the construction of P is reduced to that P + and 
P_. Each of these (say P+), can in turn be viewed as the 
composition of matrices 


p _ pO p'YTn.a.x—d^ p-fmaa: pimax-d.^ p2d'y pdj 


(A2) 


where Pp is the matrix describing how AQS dynamics 
maps the inherent structures of the landscape relative to 
7 = 7 * — dy into the set of the structures of the landscape 
associated to 7 = 7 *. The arrows in the subscript indi¬ 
cate whether P 1 2 3 * * * 7 is associated to an increase of 7 (—>•) 
or a decrease (<—) in strain. For instance, P 7 . describes 
how AQS maps the structures of the 7 = 7 * + dy land¬ 
scape onto those associated to 7 = 7 *. As the landscape 
is assumed to have M inherent structures no matter the 
value of 7 , each of these P 7 . , P 7 . is a square matrix. 
To construct each of the P2> one uses the assumption 
that the probability per unit strain to destabilize a given 
inherent structure is equal to r. So, when the strain is 
incremented by dy, a system in a given inherent struc¬ 
ture has probability 1 — rdy to be in a structure i that is 
not destabilized, and thus it maps to the same structure 
Ri in the deformed landscape through P 7 , . In that case 
the matrix element P\ ti = 1. The system has also a 


probability rdy to be in a structure Rj that is destabi¬ 
lized by the strain increment, so that it falls onto some 
randomly chosen inherent structure Rk of the deformed 
landscape, in that case the matrix element Pj) fc7 - = 1. 
Incidentally, for each configuration R 4 that is destabilized 
at 7 * as strain is increased, another configuration R m is 
correspondingly created at that strain value. This means 
that R m will be destroyed at 7 * when the strain will be 
decreased (due to the symmetry of the landscapes). This 
is a constraint on the form of the matrix P 7 . : P 7 . must 
be constructed by taking into account that the structures 
that are destroyed at 7 * when incrementing 7 are exactly 
those created at 7 * when decrementing 7 . 

The procedure outlined above is used to create a matrix 
P+ associated to some y m ax■ The P' + corresponding to 
another y' max is simply given by 


' pO p2ma.n;—dl pl ma ,n plma^-dl p2d~f pd'y 


^ 7 max < Imax 


K = 


pO P'ymasc-d-l P’lman: pima.x+d~1 pfflnax plrn.a,x plmax 

pimax+d’y pj-max pjmax—d'y p2d-y pd'y 


if irnax > 7r> 


Using this procedure its thus possible to create a se¬ 
quence of matrices relative to multiple values of 7 max- A 
similar procedure can be followed for P_). 
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FIG. 7. Single memory: (Top) KA model as discussed 
in Ref. [351 (Center) NK model as discussed in Ref. 1351 
(bottom) Fraction of inherent states that are not invariant 
under the application of a P-y r , starting from a pool of states 
trained by a different number of applications of the matrix P 71 
with 71 = 60, as a function of -y r . Data are obtained within 
the TM model setting M = 10000. It’s clear how trainings 
of increasing length yield samples that show a memory of the 
training amplitude. 
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FIG. 8. Strain profile applied to a sample to encode a double 
memory. 
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FIG. 9. (Top) KA model: Mean squared displacement be¬ 
tween configurations before and after a full deformation cycle 
of amplitude -y r , for a different number of training cycles al¬ 
ternating training amplitudes 71 = 0.06 and 71 = 0.04, as a 
function of ^ r . Data are relative to initially undeformed KA 
samples with N = 4000 and whose effective temperature is 
T = 0.466 (in reduced units see [3S] for further details). In 
the case of the longest trainings samples show a memory of 
both the training amplitudes. (Center)ATf model: Hamming 
distance between between configurations before and after a 
full deformation cycle of amplitude 7 ,., for a different number 
of training cycles alternating training amplitudes 71 = 0.3 and 
71 = 0.2, as a function of 7 ,.. Data are relative to initially 
undeformed NK samples with N = 20 and whose effective 
temperature is T = 1.0. In the case of the longest train¬ 
ings samples show a memory of both the training amplitudes. 
(Bottom) TM model: Fraction of inherent states that are not 
invariant under the application of a P-tr, starting from a pool 
of states trained by a different number of applications of the 
matrix P 71 with 71 = 60 and of P 72 with 72 = 40, as a func¬ 
tion of 7 ,.. Data are obtained within the TM model setting 
M = 10000. It’s clear how trainings of increasing N cyc pro¬ 
duce samples that show a memory of the training amplitudes. 




























